home *** CD-ROM | disk | FTP | other *** search
/ WINMX Assorted Textfiles / Ebooks.tar / Text - Mathematics - Numerical Mathematics and Computing (F).zip / amrk.f next >
Text File  |  2002-06-11  |  3KB  |  122 lines
  1. C
  2. C PAGE 406-408: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
  3. C
  4. C FILE: AMRK.FOR
  5. C
  6. C ADAMS-MOULTON METHOD FOR SYSTEMS OF ODE'S (AMRK,RKSYS,AMSYS,XPSYS)
  7. C
  8.       DIMENSION X(6)
  9.       DATA T/1.0/, X/1.0,2.0,4.0,2.0,7.0,6.0/   
  10.       DATA N/6/, H/0.01/, NSTEP/10/   
  11.       CALL AMRK(N,T,X,H,NSTEP)
  12.       STOP
  13.       END 
  14.   
  15.       SUBROUTINE XPSYS(X,F) 
  16.       DIMENSION X(6),F(6)   
  17.       F(1) = 1.0  
  18.       F(2) = X(3) 
  19.       F(3) = X(2) - X(4) - 9.0*X(3)**2 + X(5)**3 + 6.0*X(6) + 2.0*X(1)
  20.       F(4) = X(5) 
  21.       F(5) = X(6) 
  22.       F(6) = X(6) - X(3) + EXP(X(2)) - X(1) 
  23.       RETURN      
  24.       END 
  25.   
  26.       SUBROUTINE AMRK(N,T,X,H,NSTEP)  
  27.       DIMENSION  X(10),F(10,5),Z(10,5)
  28.       M = 1       
  29.       PRINT 6,T,H 
  30.       PRINT 7,(X(I),I=1,N)  
  31.       DO 2 I=1,N  
  32.         Z(I,M) = X(I)       
  33.    2  CONTINUE    
  34.       DO 3 K = 1,3
  35.         CALL RKSYS(N,M,T,Z,H,F)       
  36.         PRINT 6,T,H 
  37.         PRINT 7,(Z(I,M),I = 1,N)      
  38.    3  CONTINUE
  39.       DO 4 K = 4,NSTEP      
  40.         CALL AMSYS(N,M,T,Z,H,F,EST)   
  41.         PRINT 6,T,H,EST     
  42.         PRINT 7,(Z(I,M),I = 1,N)      
  43.    4  CONTINUE
  44.       DO 5 I=1,N  
  45.         X(I) = Z(I,M)       
  46.    5  CONTINUE
  47.       RETURN      
  48.    6  FORMAT(2X,'T,H,EST:',3(1X,E10.3)) 
  49.    7  FORMAT(8X,'X:',5(5X,E22.14))
  50.       END 
  51.   
  52.       SUBROUTINE RKSYS(N,M,T,X,H,F)   
  53.       DIMENSION  X(10,5),XP(10),F(10,5),F2(10),F3(10),F4(10)
  54.       MP1 = 1 + MOD(M,5)    
  55.       H2 = 0.5*H  
  56.       CALL XPSYS(X(1,M),F(1,M))       
  57.       DO 2 I = 1,N
  58.         XP(I) = X(I,M) + H2*F(I,M)    
  59.    2  CONTINUE
  60.       CALL XPSYS(XP,F2)     
  61.       DO 3 I = 1,N
  62.         XP(I) = X(I,M) + H2*F2(I)     
  63.    3  CONTINUE
  64.       CALL XPSYS(XP,F3)     
  65.       DO 4 I = 1,N
  66.         XP(I) = X(I,M) + H*F3(I)      
  67.    4  CONTINUE
  68.       CALL XPSYS(XP,F4)     
  69.       DO 5 I = 1,N
  70.         X(I,MP1) = X(I,M) + H*(F(I,M) + 2.0*(F2(I) + F3(I)) + F4(I))/6.0      
  71.   5   CONTINUE
  72.       T = T + H   
  73.       M = MP1     
  74.       RETURN      
  75.       END 
  76.   
  77.       SUBROUTINE AMSYS(N,M,T,X,H,F,EST) 
  78.       DIMENSION  X(10,5),XP(10),F(10,5),SUM(10),CP(4),CC(4) 
  79.       DATA  CP/55.0,-59.0,37.0,-9.0/  
  80.       DATA  CC/ 9.0, 19.0,-5.0, 1.0/  
  81.       MP1 = 1 + MOD(M,5)    
  82.       CALL XPSYS(X(1,M),F(1,M))       
  83.       DO 2 I = 1,N
  84.         SUM(I) = 0.0
  85.    2  CONTINUE
  86.       DO 4 K = 1,4
  87.         J = 1 + MOD(M-K+5,5)
  88.         DO 3 I = 1,N
  89.           SUM(I) = SUM(I) + CP(K)*F(I,J)
  90.    3    CONTINUE  
  91.    4  CONTINUE
  92.       DO 5 I = 1,N
  93.         XP(I) = X(I,M) + H*SUM(I)/24.0
  94.    5  CONTINUE    
  95.       CALL XPSYS(XP,F(1,MP1)) 
  96.       DO 6 I = 1,N
  97.         SUM(I) = 0.0
  98.    6  CONTINUE    
  99.       DO 8 K = 1,4
  100.         J = 1 + MOD(MP1-K+5,5)
  101.         DO 7 I = 1,N
  102.           SUM(I) = SUM(I) + CC(K)*F(I,J)
  103.    7    CONTINUE  
  104.    8  CONTINUE
  105.       DO 9 I = 1,N
  106.         X(I,MP1) = X(I,M) + H*SUM(I)/24.0       
  107.    9  CONTINUE    
  108.       T = T + H   
  109.       M = MP1     
  110.       DLTMAX = 0.0
  111.       DO 10 I = 1,N 
  112.         DLT = ABS( X(I,M) - XP(I) )   
  113.         IF(DLT .GT. DLTMAX) THEN      
  114.           DLTMAX = DLT      
  115.           IDLT = I
  116.         END IF    
  117.   10  CONTINUE    
  118.       EST = 19.0*DLTMAX/( 270.0*ABS(X(IDLT,M)) )
  119.       RETURN      
  120.       END 
  121.  
  122.